### mortalidade de homens 80-84 em 10 paises europeus: amostras e
### estimativas Bayes Hierárquicas

graphics.off()
windows(record=TRUE)

set.seed(180414)

library(dplyr)

## taxas verdadeiras de mort. para homens 80 a 84 em 2010, 
##   Polônia, Hungria, Reino Unido, Suécia, Rússia, 
##   Chile, EUA, Estônia, Itália, Japão
## Baixadas de http://mortality.org 

taxas = data.frame( 
  pais  = c('POL', 'HUN', 'UK',  'SWE', 'RUS', 'CHL', 'USA', 'EST', 'ITA', 'JAP'),
  mu    = c(.1064, .1157, .0772, .0779, .1319, .0914, .0732, .1073, .0784, .0712)
)

## usando dplyr (função mutate) para gerar resultados aleatórios e calcular
##   as estimativas locais de MV
X = mutate( taxas,  
               n=sample(c(50,100),10,replace=TRUE), 
               o=rpois( 10, n*mu), 
               est.local = o/n) %>% 
            arrange(mu)

X
##############################################

## hiper-a prioris
min.mu  = .05
max.mu  = .15
min.sig = 0
max.sig = .05


## log a posteriori
logP = function(mu, sig, theta) {
  ruim = (mu < min.mu) | (mu > max.mu) | (sig < min.sig) | (sig > max.sig) | (any(theta < 0))
  if (ruim) {
      return(-Inf)
  } else {
    tau  = (1/sig)^2 
    logL = -X$n * theta + X$o * log(theta)   # veross. Poisson
    logf = -1/2 * tau * (theta - mu)^2       # a priori normal
    return( sum(logL) + sum(logf))  
  }
} # logP

###############################################################
# amostragem Metropolis num espaço de 12 parâmetros 
# (mu,sigma, theta1...theta10)
###############################################################

# número de simulações
nsim = 100000 

amostra = matrix(0, nsim, 12, 
                 dimnames=list(NULL, c('mu','sig',paste0('theta',1:10))))

amostra[1,]   = c( runif(1,min.mu,max.mu), runif(1,min.sig,max.sig), rep(.10,10))
logP.atual     = logP(amostra[1,'mu'], amostra[1,'sig'], amostra[1,3:12])

aceito = 0

for (i in 2:nrow(amostra)) {
  amostra.propo = amostra[i-1,] + runif(12,-.02,.02) 
  logP.propo  = logP(amostra.propo['mu'], amostra.propo['sig'], amostra.propo[3:12])
  logu        = log(runif(1,0,1))
  if (logu < logP.propo - logP.atual) {
    amostra[i,]   = amostra.propo
    logP.atual = logP.propo 
    aceito     = aceito + 1
  }
  else {amostra[i,] = amostra[i-1,]}
}

## mantem cada 5a observ após a 1000a
mantem = seq(1000, nrow(amostra), 5)
amostra = amostra[mantem,]

##################################################
#   gráficos e resultados
##################################################

for (j in colnames(amostra)) { 
  plot( amostra[,j], main=paste('Realizações de',j), type='s',ylim=c(0,.3),col='orange')
}

for (j in colnames(amostra)) { 
  plot( density(amostra[,j]), xlim=c(0,.30),
        main=paste('Dist. a posteriori de',j))
}

resumo = apply( amostra, 2, summary)

mediana.a.post = apply(amostra,2,median)

plot(seq(-.05,.25,length=10), 1:10, type='n', xlab='theta', ylab='',bty="n",
     axes=FALSE)
axis(1)
abline(v=0)
abline(h=1:10,lty=2,col='grey')
text(-0.05,1:10, paste0(X$pais,'(',X$n,')'), adj=0, cex=.80)
points(mediana.a.post[3:12], 1:10, pch='|', cex=1.3, col='red')

conf90 = apply(amostra[,3:12], 2, quantile, c(.05,.95))
segments( conf90['5%',], 1:10, conf90['95%',], 1:10, lwd=2, col='red')

points(X$mu, 1:10, pch=16, cex=1.5)
points(X$est.local, 1:10, pch=1, cex=1.3)

list(
  erro.medio.MV    = sqrt( mean((X$est.local - X$mu)^2)),
  erro.medio.Bayes = sqrt(mean((mediana.a.post[3:12] - X$mu)^2))
)